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MicroRNAs are small noncoding RNAs that regulate genes post- 
transciptionally by binding and degrading target eukaryotic mR- 
NAs. We use a quantitative model to study gene regulation by 
microRNAs and compare it to gene regulation by prokaryotic 
small non-coding RNAs (sRNAs). Our model uses a combina- 
tion of analytic techniques as well as computational simulations 
to calculate the mean-expression and noise profiles of genes regu- 
lated by both sRNAs and microRNAs. We find that despite very 
diff"erent molecular machinery and modes of action (catalytic vs 
stoichiometric), the mean expression levels and noise profiles of 
microRNA-regulated genes are almost identical to genes regu- 
lated by prokaryotic sRNAs. MicroRNAs suppress noise when 
proteins are expressed at low levels but substantially increase 
noise at intermediate and high expression levels. This suggests 
that microRNAs and sRNAs may represent an example of con- 
vergent evolution. We extend our model to study crosstalk be- 
tween multiple mRNAs that are regulated by a single microRNA 
and show that noise is a sensitive measure of microRNA-mediated 
interaction between mRNAs. This suggests a new experimental 
strategy for uncovering the microRNA-mRNA interactions and 
testing the competing endogenous RNA (ceRNA) hypothesis. 

Author Summary 

MicroRNAs and small non-coding RNAs (sRNAs) are small RNA molecules that post-transcriptionally 
regulate target genes by binding to messenger RNAs (mRNA) and preventing translation. MicroRNAs 
are found in most eukaryotes and interact with mRNA catalytically (ie a single microRNA molecule 
regulates multiple mRNA molecules). In contrast, sRNAs are usually found in prokaryotes and regulate 
target mRNAs in a stoichiometric (one-to-one) fashion. We compare and contrast gene regulation by 
sRNAs and microRNAs using a mix of computational and analytic techniques. Despite the aforemen- 
tioned differences, we find that both mean expression levels as well as the intrinsic noise profiles of genes 
regulated by microRNAs and sRNAs are almost identical. Thus, microRNAs and sRNAs may repre- 
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sent an example of convergent evolution. We show that gene expression noise is extremely sensitive to 
microRNA-mediated interactions between niRNAs. This suggests a new experimental strategy for testing 
the competing endogenous RNA (ceRNA) hypothesis which posits that microRNAs play a key role in 
mediating interactions between mRNAs in eukaryotes. 

Introduction 

MicroRNAs are short sequences of RNA (~22 base pairs) that post-transcriptionally regulate gene ex- 
pression in eukaryotes by destabilizing target mRNAs (l][2]. Since their discovery almost two decades 
ago [3], there has been a steady increase in the number of discovered microRNAs. MicroRNAs play an 
important role in many biological processes, including animal development (4[|5], tumor suppressi on [6{[7 , 



synaptic development [8l^, programmed cell death 10 11 , and hematopoietic cell fate decisions |12l|13 



In prokaryotes, an analogous role is played by an important class of small non-coding RNAs (antisense 
sRNAs) that also act post-transcriptionally to negatively regulate proteins. These 100 bp antisense sR- 
NAs prevent translation by binding to the target mRNAs in a process mediated by the RNA chaperone 

While both microRNAs and sRNAs play similar functional roles, they act by very different mechanisms 
[17|. Eukaryotic MicroRNAs undergo extensive post-processing and are eventually incorporated into the 



RISC assembly 18 19 . The RISC complex containing the activated microRNA binds mRNAs and 
targets them for degradation. A single RISC complex molecule can degrade multiple mRNA molecules 
suggesting that microRNAs act catalytically to repress protein production. In contrast, both the mRNAs 



and sRNAs are degraded when bound to Hfq 15 16 suggesting that prokaryotic sRNAs, unlike their 



eukaryotic counterparts, act stoichiometrically on their targets. 

While stoichiometric regulation has been extensively studied theoretically and experimentally [20 
there exist relatively little work on understanding microRNAs and other catalytic non-coding RNAs 
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29 . Both theoretical calculations (see Supporting Information of [20]) and quantitative experiments 27 
indicate that mean protein-expression of microRNA regulated genes follows a linear-threshold behavior 
similar to that of sRNAs in prokaryotes. In contrast, it was argued in ^29^ that the intrinsic noise 
properties of catalytic sRNAs/microRNAs differ significantly from sRNAs. Presently, it is unclear how 
to reconcile these results and it raises the natural question of how the differing mechanisms employed by 
sRNAs and microRNAs affect the intrinsic noise profiles of regulated genes. 

To answer this question, we used a generalized model of gene regulation by non-coding RNAs to 
calculate the mean expression and intrinsic noise of regulated proteins. Our model is based on stochastic 
mass-action kinetics with tunable parameters that allow us to vary the mode of action from stoichiometric 
interactions such as those in prokaryotic sRNAs to highly catalytic interactions that mimic eukaryotic 
microRNAs. Contrary to [29| , we show that in many parameter regimes the intrinsic noise properties 
of microRNAs and sRNAs are qualitatively similar. Finally, motivated by the renewed interest in the 
competing endogenous RNA (ceRNA) hypothesis which suggests that microRNAs induce an extensive 
mRNA-interaction network, we extended our model to consider the case where a microRNA regulates 
multiple mRNAs. We calculate the intrinsic noise for these models and show that noise is an extremely 
sensitive measure of cross talk between mRNAs. This suggests a new experimental strategy for uncovering 
microRNA-induced mRNA interactions. 

Results 

Model description 



Here we propose a model of gene regulation by non-coding RNAs based on mass-action kinetics 29 . A 



schematic of the model is shown in Figure [T]\. Our model has four species: mRNA molecules denoted 
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by TO, functional non-coding RNAs denoted by s, mRNA- noncoding RNA complexes denoted by c, 
and the number of expressed proteins denoted by p. We note that s can represent either the number of 
prokaryotic small RNAs or the number of functional microRNAs found within the RISC complex. mRNA 
molecules are transcribed at the rates am and translated at a rate Up. In the absence of regulation, 
mRNAs degrade at a rate t~^. Ug represents the mean rate at which functional non-coding RNAs are 
formed. For prokaryotic sRNAs, this is simply the transcription rate of sRNAs. For microRNAs, this 
is an aggregate rate that accounts for the complicated kinetics of transcription and assembly into the 
RISC complex. mRNAs and noncoding RNAs form complexes c at a rate /i and disassociate at rate 
7. Importantly, the complexes can also be degraded at a rate t^^. This degradation can be actively 
regulated or conversely stem from dilution due to cell division. 

Once mRNAs bind noncoding RNA and form a complex c, they can no longer be translated, resulting 
in decreased protein expression. To account for the differences between microRNAs and sRNAs, we have 
an additional parameter /3 which measures the "recycling rate" of noncoding RNA. When /3 ^ r^^, 7, the 
non-coding RNAs function catalytically with a single non-coding RNA molecule able to bind and degrade 
multiple mRNA molecules. This is a good model for gene regulation by microRNA in eukaryotes. In 
contrast, when /3 ^ the noncoding RNAs function stoichiometrically. In particular, for prokaryotic 
sRNAs it is commonly believed that /3 = and no recycling of noncoding RNAs takes place. Finally, we 
note that the ratio of /3 to 7 is a measure of how much of the gene regulation happens through nucleolytic 
cleavage as opposed to translational repression. 

We use two different approaches to mathematically model gene regulation by non-coding RNAs. 
We calculate both the mean expression levels as well as "intrinsic noise" due to stochasticity in the 
underlying biochemical reactions. First, we perform simple Monte-Carlo simulations for the reaction 



scheme outline above using the Gillespie algorithm 30 . While these simulations are exact, this method 



is computationally intensive making it difficult to systematically explore how noise properties depend 
on kinetic parameters. For this reason, we use a second approach based on linear noise approximation 



(LNA) [22 31 32 . The LNA approximates the exact master equations using Langevin equations that 
correctly reproduce means and variances. This allows us to derive analytical expressions for the noise 
and systematically explore how gene expression noise depends on parameters. 

Within the LNA, we can mathematically represent our model using the equations 

ds _i „ 

— = as-T^ s + pc + ^c- fims + r]s + rip + rj^ - 77^ 
dm _i 

— = am- T„ TO + 7c - fims + rim + ri^ - Vt^ 

<" 

— = fims - pc- ^c- c + rjc-rip-ri^ + rj^ 
dp _i 

= rv^m. — T 

dt 



ttpTO - Tp p + rip 



with rjsjTjmT'nci and rjp, being the noise in the the birth-death processes of non-coding RNAs, mRNAs, 
complex, and proteins respectively, 77^ the binding noise, 77^ the unbinding noise, and r]p the RNA recycling 
noise. The variance of these noise terms is given by 



('7i(^))=0, i = s,l3,-/,fi,m,c,p 
\j(t')) = 5ijC 



{m{t)i^,{t'))^6,,a}5{t~t') 



(2) 



where 5ij is the Kronecker delta and 



2 I —1- 2 0- 2 - 2 /o\ 

a^^ as+ Tg s, o-o = Pc, a = 7c, a fims, (3) 



<^m ^ am + Tm^m, ^ T^'^C, = apTTi + Tp'^ p, 
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with n the steady-state concentration of species n. We modeled each interaction in Figure [T|\ as an 
independent Poisson process with a mean rate given by mass-action kinetics fsTl . Each Langevin term is 
an independent process. Care must be taken to ensure the correct sign in the cross-correlations [32|. In 
the remainder of the paper, analytical results from the linear noise approximation are shown along with 
simulations using the Gillespie algorithm [30' . Both methods are in good agreement. 

Mean Expression Levels 

We begin by deriving the steady-state response of our system by setting the time-derivatives of the left 
hand side of equation ([l]) to zero. Denoting the steady-state concentration of a species n by h, we find 
that 

(/jofj, - a„i - \ + J {(f)as - a,n - A)^ + 4X(j>as 
s — 1 

Oim - (t^ois - >^ + y {am - 4'as - Xf + 4Xa„i (4) 
m = — 

C — flfhSTcR 

p = apTpfh 
where we have defined the quantities 

'/'=l+/3r,, (5) 

A=^^, (6) 



and 



r-^^=/3 + 7 + T--i (7) 



'cR 



P + J + Tc 



(8) 



For all values of the parameters, the steady-state protein levels exhibit a linear-threshold behavior in 
mean protein concentration (Figure IB) with the threshold set by the condition am ~ (jjcfs- Such linear- 
threshold behavior in mean protein production has been extensively studied in the context of non-coding 
RNAs and it is useful to divide gene expression into three distinct regimes 20 22 27 . In the repressed 
regime (a™ ^ ^ctg) non-coding RNAs almost always bind mRNAs and prevent translation. In contrast, 
in the expressing regime 3> (j^a^) there are many more mRNAs then noncoding-RNAs resulting in 
protein production that varies linearly with am ■ Finally, there is a crossover between these two behaviors 
when am ~ <pas- 

The factor <p multiplying as renormalizes the transcription rate of the non-coding RNA to account 
for the fact that single non-coding RNA can degrade multiple mRNAs. To see this note that 4>~^ = 
T~^/(/3 + T^^) is just the probability that a non-coding RNA is degraded when it is bound to an mRNA 
in a complex under the assumption complexes do not disassociate( 7 <C /3 or 7 ^ t~^). Thus, we 
can think of (p as the average number of mRNAs that a non-coding RNA will degrade before it is itself 
degraded. As expected, when /? = 0, = 1 and these results reduce to those derived for prokaryotic 
sRNAs 20 22 . Overall, this linear threshold behavior with variable rescaling is consistent with recent 



experiments on microRNAs 27 
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Intrinsic Noise 

Gene regulation is inherently stochastic due to the small number of molecules present in the cell. Noise 
has important phenotypic consequences for a variety of biological phenomena '33' and for this reason it 
is important to characterize the intrinsic noise properties of non-coding RNAs. As is usual, we define 

the intrinsic noise as the variance in protein level divided by mean protein level squared, ^''^ ^^^^P ^ ; 
where brackets represent steady state time averaging. This is a measure of relative protein fiuctuations 
compared to their mean. Study of noise in bacterial sRNAs shows a peak in the crossover regime that 
emerges as a result of competition between mRNA and sRNA [22). The switch-like behavior of the 
system, due to its linear-threshold nature, results in an amplification of any fluctuations in the mRNA 
level that is in excess of the mRNAs bound to noncoding-RNA. As a result, noise is increased in the 
crossover regime. We observe a similar qualitative behavior for noise in all parameter regimes of our 
model. Figure [2] is a plot of protein noise as a function of mean protein concentration for catalytic and 
stoichiometric interactions, showing a peak at protein levels that correspond to the crossover regime in 
both cases. The height of the peak increases with the binding affinity /i of a non-coding RNA for its 
target mRNA. On a whole, the noise profile of catalytic and non-catalytic genes is remarkably similar. 
There are however slight differences. The peak is slightly higher and occurs at a slightly larger mean 
protein level for catalytic non-coding RNAs. This suggests that the underlying reason for the noise peak 
is not the enzymatic mode of action of non-coding RNAs, but the fact that the number of mRNAs and 
the number of non-coding RNAs (appropriately normalized by (/)) are almost equal. 

To better understand the effect of catalytic interaction on noise, we calculated protein noise versus 
am for different values of /3. Figure [Sj'V shows the results for Gillespie simulation and linear noise 
approximation. To meaningfully compare noise between stoichiometric (i.e. prokaryotic sRNA, /3 = 0) 
and catalytic (i.e. eukaryotic microRNA, /3 ^ T~^,'y) we need to compare noise at the same mean 
protein level. Figure [3j3 shows this comparison for noise in the repressed regime as a function of f3 with 
am chosen such as to keep mean protein concentration constant. As can be seen, the noise decreases from 
its original value in the stoichiometric regime (/? = 0) to a slightly lower value as (3 is slightly increased 
(/? « 2) and any further increase in /3 does not affect noise dramatically. We also plotted the maximum 
of the noise peak in the crossover regime using both linear noise approximation and Gillespie algorithm 
(Figure [3p) . Notice that the peak height initially increases as a function of P and then reaches a plateau 
for large /? upon entering the catalytic regime. 



Including transcriptional bursting 



Experimental evidence suggests that mRNAs are often produced in bursts 34 . Transcriptional bursting 



represents another important source of stochasticity that was ignored in the analysis presented above. 
We can extend the model presented above to incorporate transcriptional bursting by considering the case 
where the genes encoding for mRNAs can be in two distinct states: an "on" state where mRNAs can 
be transcribed and an "off" state where transcription is not possible. For example, in eukaryotes the 
two states may correspond to whether the chromatin is condensed or not 35 . To model transcriptional 
bursting we replace the equation for mRNA production in Eq. ([I]) by the pair of equations 



dg 
dt 
dm 
~dt 



^ k^{l - g) - k+g + r]g 



(9) 



+ 7C - fims + rjm + 111 ~ Vt^ 



where the probability that a gene is in the on state is denoted by g, k^{k+) are the on(off) rate of the 
gene, and is the maximum mRNA transcription rate. The mean probability of the gene being on is 
given hy g — k^/ (k^ + 22 . The gene noise rjg is Gaussian white noise with mean zero and variance 



6 



given by {rig{t)rig{t')) = 2k+g5{t-t'). The mRNA noise rj^ is now {T]m{t)Vm{t')) = ia^g+T;;^^m)S{t-t'). 
All other equations remain the same as before. 

In the analysis that follows, we assume that /c_ is fixed but fc+ can vary. This corresponds to the 
biological situation where a gene is regulated by a repressor that can turn the gene off. To compare noise 
for different values of /?, we choose so that the mean protein levels remain constant. Furthermore, we 
concentrate on the repressed regime (see Figure |4]) . Here noise decreases slightly as (3 is increased which 
is very similar to the case without bursting as was shown in Figure [3j3. This means that noise in the 
repressed regime is insensitive to bursting regardless of how catalytic the interaction is. 



Asymptotic Formulas for Noise 
Expressed 

To gain further insight, we have derived asymptotic formulas for the noise in the repressed and expressing 
regime. In the expressing regime with large mRNA transcription, we define the small parameter e = 
^ with A given by Eq. 6 In the expressing regime, e ^ 1, and the steady-state expression levels 

of mRNA and non-coding RlNTAs take the form 

m~ , s ^ agTge (10) 

qfiTs e 

with q given by Eq. |8] Furthermore, the protein noise in this regime is identical to that of a transcrip- 
tionally regulated gene and is given by 

{(spr) _ .1+6 
p^ p p 

where h = apTm (see Materials and Methods) and it is assumed that mRNAs degrade much faster than 
proteins (rp ^ t^)- The quantity 6 is often called the 'burst size' and measures the average number of 
proteins made from a single mRNA molecule [36]. This result shows that in the expressing regime the 
protein noise has no /3 dependence, and thus the protein noise of post-transcriptional and transcriptional 
regulation are identical in the asymptotic limit. 



Repressed 

In the repressed regime, we now have e = '^^ <C 1 and the average number of mRNA and non-coding 
RNA molecules is given by 

m ~ arriTme (11) 

1 



qpTmfpe 

where q 

^ /3 + 7 + T,-i' 

is the probability that a complex is degraded. The noise in this regime is given by 

(Ml ^ l±c^ (12) 

p 

where C is a constant that is dependent on both /3 and 7 (see Materials and Methods). When /3 or 
^ 7 it can be shown that C « 1. Note that this condition includes the catalytic regime with 
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/3 ^ 7i ''"tT^ ■'^'^U as the stoichiometric regime /? = 0, r^"'^ ^ 7 (see Methods). Thus, we conclude that 
in the repressed regime, non-coding RNAs reduce noise. In particular, proteins are now produced from 
mRNAs in a smaller burst with typical size given by e^b <IC b. Since the burst size is just the average 
number of proteins made from an mRNA, b = afcr™, we can equivalently interpret this results as changing 



the effective lifetime of the mRNA molecules from t„i to e(Tm 20 22 



Scaling Behavior Near Crossover Regime 

Our analysis show that the width of noise peak at the crossover regime is independent of recycling ratio 
j3. To understand this behavior better, we studied the crossover regime for different parameter values 
using linear noise approximation. Figure j5j shows the protein noise and protein mean as a function of 
after rescaling both by their value at the the crossover {am = 4"^s) for various recycling ratios. Notice 
that for 7 ^ r^^ these normalized plots of protein mean and noise show an approximate data collapse. 
As 7 is increased, this scaling behavior breaks down (Figure [5]) . The collapse of data for mean protein 
can be analytically derived given the fact that mean protein is only dependent on /3, 7, Tc through the 
combination qd) = ^^'^^ 



This scaling of the mean protein number can be better understood if we define x = and define 
p{x) as the mean number of proteins corresponding to this value. Dividing this quantity by its value at 
X = 1 and using equation ^ results in the following equation for the normalized mean protein: 



Ml) 



x-l-^ + \l{x-l-^]\4^x 



(13) 



,, , , + 4-^ 

90 V \Q<t> J 90 



where v = — is a constant with no dependence on /3,7, Tj,. Thus, the normalized mean protein 
number depends on /3,7,Tc only through the combination qcj) = ■ The parameter q(j) is the 

probability a complex will disassociate. In the limit 7 <C P,t~^ this parameter will be equal to 1 and the 
scaled mean protein level becomes independent of j3 causing the curves for different /3 to collapse on top 
of each other (figure [5]4_) . It is also interesting to note that any other condition on the parameters that 
removes the dependence of q(t> on j3 will also have the same effect (e.g. 7,/? ^ Tc"^)- Somewhat more 
surprisingly, the noise also shows an approximate data collapse. We suspect that this collapse is likely 
indicative of universality near the crossover between the repressed and expressing regimes. 

mRNA Crosstalk and the ceRNA Hypothesis 

Recently, the competing endogenous RNA (ceRNA) hypothesis proposed that microRNAs may play a 
crucial role in the cell in global gene regulation by inducing interactions between mRNA species [37] . 
The central mechanism underlying the ceRNA hypothesis is the idea that mRNA species may have 
interactions amongst themselves that are not direct but are instead indirect and mediated by competition 
and depletion of shared microRNA pools. The hypothesis is that these indirect mRNA interactions results 
in a biologically important mRNA network. However, the breadth and strength of microRNA induced 
interactions in eukaryotic genomes is still not well understood. For this reason, it is crucial to develop 
new strategies for measuring microRNA induced interactions between commonly regulated mRNAs. To 
this end, we asked whether the noise profile of regulated mRNAs could be used to uncover microRNA 
induced interactions. As a first step, we studied the simplified case where two different mRNA species 
are regulated by a single microRNA and compete over a common pool of microRNAs and focused on the 
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effect of microRNA-induced crosstalk between mRNA species on the noise properties of regulated genes. 
The results presented here can be easily generalized to the case of many mRNAs interacting with many 
microRNAs. 

A schematic of our simplified model is shown in Figure [6j^. Two species of mRNAs are regulated by 
a common microRNA. Notice that the mRNAs do not directly interact in the model and all interactions 
are indirectly induced by the shared microRNA pool. We can represent this using a straight-forward 
generalization of the model considered earlier 

Cts - T^^S + PiCi + /32C2 + 7lCl + 72C2 - (MI?"! + M2™2)s + ?7s + Tjp^ + Tjp^ + 7/^^ + 7/^^ ^ Vl-ii ^ 

(14) 

^J,im^s - PiCi - y,Ci - T~\i + rid - ?7/j, - r/^, + 7?^^ 
apWii - Tp^p + rip 
with 

{r]j{t))=0, j = s,(3i,ji,fj,i,mi,Ci,p ,« = 1,2 (15) 

and variances reflecting the Poisson nature of interactions: 

= + T~^s, afj. = P.Ci, a\ = jiCi, cr^^ = f^im^s (16) 
CT,^„^ = am, + T;^]mi, al. = T^^^Ci, = Upfhi + r^^p , j = 1, 2 

The binding of microRNAs to mRNAs is controlled by the interaction rates /ii,2- These rates reflect 
the binding affinity of microRNAs for the two mRNA species. We asked how protein noise and means 
change as we vary the transcription rates, am-i a: of the mRNAs. The remaining parameters are assumed 
to be the same for both mRNAs and we have suppressed the indices on these parameters. MicroRNAs 
function catalytically and we focus on the parameter regime /3 ^ T'^^,"f. Figure [6|3 and C shows the 
mean protein levels, pi, and intrinsic noise of protein species 1 as a function of the the transcription rates 
of the two mRNA genes, ami 2j the case of equal binding affinities (/ii — /i2). Notice there is a peak 
in the noise similar to the single-species case. Once again there is good agreement between simulation 
and analytic calculations based on Langevin noise. We also examined the case where the mRNAs have 
different binding affinities for the microRNA, /zi = 0.2, ^2 = 2. This results in an asymmetry in the noise 
peak but does not change the major qualitative features of our results (see Figure [6|l) and E) 

As in the single-mRNA species case, the behavior of the system can be divided into regimes depending 
on whether the combined transcription rate of both mRNA species is bigger or smaller than the normalized 
sRNA transcription rates. The crossover regime, + ctm^ ~ 't"^s splits the transcription rate space 
(amijama) i^i^o an expressing regime, -I- > 'P^'-s, and a repressing regime, + 0,^2 < 4"^s- 
Here, similar to the single species case we see a sharp peak in the noise at the crossover regime. Since 
the cross-over regime depends on the total transcription rate, ami +am2, this suggests a general strategy 
for uncovering microRNA interactions based on "intrinsic noise spectroscopy". To see if an mRNA, 
say species 2, interacts with an mRNA species 1 that is repressed by a microRNA {ami ^ 4"^s), we 
can monitor the intrinsic noise of protein 1 (i.e. a peak) as we vary am2 while keeping ami fixed. If 
mRNA species 2 interacts with the microRNA, there will be a dramatic signature of the interaction in 
the intrinsic noise of protein 1. In contrast, there are no dramatic features in the mean levels of proteins. 
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Discussion 

In this work, we studied gene regulation by non-coding RNAs. Whereas gene regulation by prokaryotic 



sRNAs has been extensively studied 20 -26 , there exist relatively few models of gene regulation by 
catalytic microRNAs 27 -29 . Here, we used a simple kinetic model to study both the mean expression 
levels and intrinsic noise properties of post-transcriptional regulation by non-coding RNAs. Using a single 
parameter, our model interpolates between the stoichiometric behavior of prokaryotic sRNAs and the 
catalytic behavior characteristic of eukaryotic microRNAs. We found that both sRNAs and microRNAs 
exhibit a linear threshold behavior, with expressing and repressed regimes separated by a crossover 
regime. At the crossover, the mRNA transcription rate roughly equals the product of the non-coding- 
RNA transcription rate and the average number of mRNA molecules degraded by a single non-coding 
RNA molecule. In all cases, the intrinsic noise was smaller in the repressed regime and showed a sharp 
peak in the crossover regime. We found that for most parameter regimes, the intrinsic noise in the 
crossover regime shows an approximate data collapse, giving hints that there may be universal behavior 
near the transition from the repressed to expressing regime. We then generalized our calculations to study 
between two mRNAs regulated by a single microRNA. We found that the intrinsic noise is an extremely 
sensitive measure of microRNA induced interactions between mRNAs. 



Our results for the mean expression profile is consistent with recent experimental studies 27 . How- 



ever, our conclusions about the intrinsic noise profiles of catalytic non-coding RNAs are different from 



Hao et al. 29 . They concluded that the intrinsic noise profiles of catalytic and stochiometric differed 
significantly. The reason for this discrepancy is that Hao et al. did not normalize the sRNA transcription 
rate as by the recycling rate 0. Consequently, they compared the cross-over regions of sRNA regulated 
genes to repressed regions of microRNAs. As shown above, after making this normalization there is 
extensive data collapse of intrinsic noise profiles for both stoichiometric and catalytic genes. 

One of the striking results of our calculation is the similarity between sRNA-regulated and microRNA- 
regulated genes. This similiarity is somewhat surprising given that microRNAs and sRNAs are found 
in different kingdoms (prokaryotes versus eukaryotes) and utilize distinct biochemistry and molecular 
machinery. Eukaryotic microRNAs use complicated nuclear machinery to export microRNA into cy- 
toplasm and bring it to mature state by incorporating the RNA strand into the RISC complex. In 
contrast, prokaryotic sRNAs undergo relatively little post-processing and bind mRNAs via the chaperone 
protein Hfq. Given these extensive differences, the similarity between the expression characteristics of 
microRNA-regulated and sRNA regulated genes are suggestive of convergent evolution. 

Our calculations show that mRNAs regulated by catalytic non-coding RNAs have large peaks in 
the intrinsic noise. This differs significantly from results that would be derived from more traditional 
treatments of catalytic interactions based on the Michaelis-Menten or Hill equations. The underlying 
reason for this difference is that traditionally, the Michaelis-Menten equations are derived under the 
assumption that substrates of enzymes are in excess compared to the enzymes themselves. In contrast, 
here we are interested in the case where the number of enzymes (microRNAs) and number of substrates 
(mRNAs) are comparable. This accounts for the sharp peak in noise observed in the crossover regime in 
our models that is absent in traditional treatments of enzyme kinetics. 



Our calculations also suggest a strategy for testing ceRNA hypothesis 37 , which posits that microR- 
NAs induce extensive interactions between mRNA molecules. Our calculations suggest that protein noise 
is a more sensitive measure of the competition between the two species than mean levels. Thus, it may 
be easier to uncover interactions between mRNA by measuring changes in the noise of regulated genes. 
We suspect that this will be true even when we generalize our calculations to consider the case where n 
mRNA species compete over the same pool of microRNAs. In this case, we hypothesize that there would 
still be a sharp peak in the intrinsic noise when the total transcription rate of all mRNAs equals the 
appropriately normalized transcription rate of microRNAs. In the future, it will be interesting to analyze 
these more complicated models in greater detail. 
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Materials and Methods 

Single Species Linear noise approximation 

Linearization of equation [T] results in: 



di 



6s' 




5m 




6c_ 





'sR 
— flffl 

fifh 



—^s /3 + 7 



^cR 





5s' 








dm 


+ 


Vrn + Vl - 7?M 




Sc_ 







where we have named the transfer matrix A and we define the renormalized hfetimes: 

-1 



—1 —1 - —1 —1 



TcM = '^c^ + P + 1 



(17) 



To find the solution of this equation we apply a Fourier transform and solve the resulting equation: 



Ss' 








-1 


6m 




^m 


+ T^R -7 




Sc 











where tilde denotes Fourier transform. So we have: 

1 



5m 



[ F{io) Giuj) H{uj)] 



{iuj — Ai)(iw — A2)(ia; — A3) 

where Ai, A2, A3 are the eigenvalues of A and: 
F{u}) — jiirh — fim{iuj + t^j^) 

H{lu) = -(/? + "f)txrfi + -f{iu + T^j^) 

F(w)r), + [F{lu) - H{io)]fip + [F{u) + G{u) - H{u)](f,^ - ry^) + G{uj)fim + i?(w)?)c 



(18) 
(19) 
(20) 



or 



[iuj — \i)[iLj — \2){iuJ — A3) 



Since 5n = '?p "p ™ have 



(21) 



Q 



\F{Lo){'ai + iFiu) - H{u)\-'ai + |F(..) + G(a.) - il(c.)|^(a^ + a^) + |G(^)|V,^ + \Hi^)\^a^ 



(a.2 + Af)(a.2 + A^)(a.2 + A^) 



We can rewrite Q as: 
Q = Xw" + Yuj^ + Z 



(22) 



with 



Y = {iimf{as + T~^s) + {iifh + ■y)'^ PrcRUfhs 
+ (7-7^ + + Pfi^rcR + l)firns 

+ [Kr - ^cr)^ + 2(/3 + 7)nrhf]iam + T-^m) + j^i^fhs 
Z = (7 - T^j^f(pfaf{as + T~'^s) + {T~'^lJ.fn + 'yT~'^f ^TcRUffis 



now if we write 

U?i\^\ = '^ 2,2 [ + Yuj"^ + Z 

\KP}) ^^p + ^pj 27r(c.2 + A2)(c.2 + A2)(c.2 + Ai)(u;2+Tp-2) 

and let A4 = —t~^ by use of partial fractions and integration we get 

i=l ' jjti j ' 

Single Species Asymptotic calculations 
Expressing Regime 

In the expressing regime with large mRNA transcription rate we demand e 
m and s in terms of e we will get: 

m = — ^ + 0(e) 
s = agTsC + O(e^) 



To find the eigenvalues of A we will expand it in terms of e: 



-rr' - rr^Q-^e-^ + 0(e) -l^a^ne + 0{e^) 



T~^q-h-^ + 0{e) 



0{e) 



the eigenvalues of this matrix satisfy the following equation: 
+ a2X^ + aiX + ao = 



where 



02 = ^q ^ 



0(1) 



a,=T-\-\T-'+T-')e-'+0{l) 
ao = g- V- V-^e-i + 0(1) 



' + 7 

7 
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this equation can be analytically solved by expanding A's in terms of e, and results in the following: 

Ai = -r„i+0(e) (30) 



A, 



0(e) 



A3 = -^^+0(l) 

which shows one fast mode and two slow modes. Next we will calculate the noise by expanding equation 
24 in terms of e, which results in: 



X 



Y 



2A 



(31) 



2A 



2A 



plugging this result into equation 25 and using the approximation Tp ^ r,„ we will have: 

(jSp)') .1 + 6 
p 

Repressing Regime 

In this regime with mRNA transcription very small, we demand e = ^ 1. Expansion of to and 

s in terms of e gives: 



TO = a„ir„ie + O(e^) s 



1 



qfj,Tjn(l)e 

Next we expand A in terms of e: 



0(e) 



A 



-Aia™T„e + 0(e2) -^--i - -i^ + 0(e) 7 



0(e) 



after some calculation we find the following closed form for eigenvalues of A: 
A^ + a2\^ + aiA + ao = 
where 

a2-r„ig-ir'e"' + 0(l) 



ai=T^ q 

ao = T,-'q-'r'r-\T-' + f3)e-' + 0(1) 



(32) 



(33) 
(34) 



this equation can be analytically solved by expanding A's in terms of e, and results in the following: 

Ai = -r.-i + 0(e) (35) 
A2--(r-i+/3) + 0(e) 



A, = - 



Oil) 
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which has one fast mode and two slow modes. Next we calculate the noise by expanding equation |24| in 
terms of e: 



-1^-1 



(36) 



+ + tJ) 



plugging this result into equation [25] gives: 
{{5pf) 1 be q^<f>^ ( 



+ 



Xt-'^ - Yt-' + Z 



X{T-^+P)^-Y{T~^+Pf + Z 
[rp^ - {tc' + /3)2) (r.-' - (rc-i + /3)2) {r,-' + /3) ^ 

for /3 7, T~^, t7^, T~^, after some calculation we get ^^"^^2^ ^ i (1 + 6e). 

ceRNA Linear Noise Approximation 

For two mRNAs, we linearized equation [T5| as: 



(37) 



dSx 
dt 





5s ' 




r --7.^ 




-M2S 


A +71 


^2 + 72 





" 




a" 




6mi 




-A*i™i 







71 

















Sm2 









'mi?,2 





72 












Sx = 


Sci 


, A = 


























SC2 









M2S 





^cR2 














Spi 






























Sp2_ 


















-<\ 






with 
























1 -1 






















+ fiiffii + 


M2TO2, 




'''mi f^iS, 




= -cT^ + 


7i + A 







and: 



where 



r = 



51 


92 


53 


54 


55 








92 


96 





57 











93 





58 





59 








94 


57 





510 











95 





59 





511 























512 























5l3 



(38) 
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and 



91 




^1 + 


4.+ 


<- 


92 




<4 


-< 




93 


= 




-< 




94 


= 




-< 




95 










96 






^92 




97 




-92 






98 






f 33 




99 




-93 






910 




<~ 


54 




911 




<- 


95 




912 










913 











'P2 



(39) 



We find corresponding two point correlation functions by use of the following equation 32 



r 



qs 



\p + Xr 



where A's are the eigenvalues of A, and B is the matrix of eigenvectors, according to: 
''^AijBij = XkBik 



we saw a good agreement between these results and Gillepie simulations. The two methods showed at 
most a deviation of 30% from each other. Due to this high similarity and to save space only the results 
of Gillespie simulations are shown in figure |6] 

ceRNA asymptotics 

Here we assume that other than /z's all the other parameters are equal between the two species. As a 
result the corresponding indices have been dropped in the following calculations. For the steady state 
values of concentrations we have: 



Us = - (/3 + 7) 
which results in 



' cR 



s + {pimi + fJ,2m2)s 



as = Ts s + q 



ami Ml 



Oim2fJ'2 



7'rn 
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After some calculations we derive the following polynomial for sRNA concentration, s: 

f + {Bi +B2 + A1+A2- K)s^ + {B1B2 - K{Bi + B2) + A1B2 + A2Bi)s - KB1B2 = 

with 

K ^ asTs, A^ = — - — , Bi ^ —— 
V <P 

furthermore in the following calculations — ^ and index T represents total sum over z, such that: 



At = Ai+A2, Bt^Bi+B2, At = Ai+A2, a^T 
Expressing regime 

In this regime At 3> G = max{BT, KBt, B1B2) and we can simplify the polynomial equation by defining 
e = G/At and multiplying it by e while keeping coefHecients to first order: 

es^ + Gs^ + GDs - KBiB2e = 

with D = + "j^Bi . Note that GD is of order ©(e") and does not require e expansion. The equation 

for s has the following asymptotic solutions: 

G ^ kBiB2 

with the only positive solution being s = ^ g^^^ £ — Qm"°^°r,i2 ■ In the limit of = 0, this reduces to the 

single species result in equation 27 Finally for mRNA we have m ~ TmOCmi(^ — I'pTmlJ'is) ^ TmCtmi +0(e) 
which is the expected result. 

Repressing regime 

In this regime K ^ G = max{Ai, Bi) which is equivalent to (t>as 3> Q;m;,Ai. Using this assumption 
we can simplify the polynomial equation by defining e = G/K and multiplying it by e while keeping 
coeffiecients to first order: 

es^ - (G - {At + BT)e)s^ + {{B1B2 + A1B2 + ^2^1)6 - GBt)s - GB1B2 = 

which has the following asymptotic solutions: 

G 

e 

with X being the solution for the following equation: 

+ xf - (G - {At + BT)e){^ + xf - GBt^ - 



This results in a; = —At- So the only positive solution is s ^ — — At ~ -it{4'Cts~oi„ij,) and 



(hOLs—Oijn^ 



which for = reduce to our single species results in equation 32 



Crossover regime 

Crossover regime is roughly where the repressing regime's line intersects with s = 0, hence at this point 
we have — (pets 
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Figure Legends 




Figure 1. microRNA model: A) Schematic of the interactions between noncoding RNAs and 
mRNA. as and a,n are respectively the transcription rate of microRNA and mRNA. r~^, r^~^, t^T^, r^"^ 
are respectively the degradation rates of mRNA, noncoding RNA, complex and protein, /i and 7 are 
respectively the direct and reverse interaction rates between mRNA and noncoding RNA, and /3 is the 
catalyticity. B) Analytical results showing protein mean versus normalized transcription rate, ^j^, for 
different values of /i in the catalytic regime (/? = 10, Tc = 1) where <j> = 1 + Ptc- The dashed line is the 
theoretical result for infinitely large /i. The following parameters have been used in this plot: 
as = 3, ttp = 4, Ts = 30, T„, = 10, Tp = 30, 7 = 1. 
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Figure 2. Protein noise in two regimes: Simulation results showing protein noise as a function of 
mean protein concentration for stoichiometric and catalytic interactions. For high /i (solid curves) there 
is higher noise in the crossover regime in the catalytic case. Same results are plotted for smaller /i 
(dashed curves) showing a decrease in noise at the crossover regime while maintaining the same overall 
behavior. Parameters same as in Figure [Tj3. Stoichiometric regime with l3 — and catalytic regime with 
/3 = 10. 
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Figure 3. EfTect of catalyticity on noise : A) Simulation results for noise as a function of /3 and 
Inset: analytical results for the same system. Parameters same as Figure m with ^ = 2. B) Noise 
in the repressed regime as a function of /3 for constant protein mean (p = 10). Inset: protein mean as a 
function of /3. C) Maximum of noise in the crossover regime as a function of /3. Inset: maximum of 
Fano factor in the crossover regime as a function of /3. 
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Figure 4. Catalyticity and Bursting: Noise in the repressed regime with bursting as a function of 
/? for constant protein mean. For each data point a^, is chosen such that < p >= 10. Furthermore 
a™ = 10, A:„ = 1, fc+ = — 1). The remaining parameters are same as in Figure jlj with /i = 2. 

Inset: protein mean as a function of /3. 
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Figure 5. Scaling at crossover regime: Parameters same as Figure [l] with fi = 2. A) Protein 
noise normalized to its value at a™ = (pa^ plotted as a function of for 7 = 0. Each line is a different 
value of /3. Same legend for all figures. B) Protein mean normalized to its value at am = 4>o^s plotted as 
a function of for 7 = 0. C,D) Graphs similar to A,B plotted for 7 = 1. E,F) Graphs similar to 
A,B plotted for '7 = 10. 




Figure 6. ceRNA hypothesis: A) Schematic of mRNA crosstalk through a shared pool of 
microRNAs. 2 stand for transcription rate of each mRNA and /ii_2 correspond to interaction rates 
between mRNA and microRNA. The other interactions are as in Figure lA. B,C) Protein mean (B) 
and noise (C) as a function of transcription rates of the two mRNAs with equal /i's. Parameters same 
as in Figure [1] with /ii = /i2 = 2. D,E) Protein mean (D) and noise (E) as a function of transcription 
rates of the two mRNAs with unequal /i's. Parameters same as in Figure [TJwith = 0.2,^2 = 2. Noise 
surface plots have been smoothed and interpolated for better visibillity 



